%% Optimize JRD model
clc
clear
rng('shuffle')

global newTestProbe
global obsRespErr
global refDirectionFull
global layout

% Prepare data
% Individual fitting
Exp2matrix = readmatrix('KellyMcNamara2010_Exp1.xlsx');

% Initial parameter values
x0 = [.4 4];

Exp2matrixL1 = Exp2matrix(find(Exp2matrix(:,4) == 1),:);
Exp2matrixL2 = Exp2matrix(find(Exp2matrix(:,4) == 2),:);

% Need to add one to each object to shift indexes
Exp2matrixL1 = Exp2matrixL1 + 1;
Exp2matrixL2 = Exp2matrixL2 + 1;

% Use this for global search
opts = optimset('Algorithm','interior-point');
problem = createOptimProblem('fmincon','options',opts,'objective',@JRDLogGrowthModelFullRefFixL_Study2, ... 
    'x0',x0,'lb',[0,0]);
gs = GlobalSearch;

for p=1:1
layout = 0;
refDirectionFull = 1;

for i=1:16
    if layout == 1
        newTestProbe = Exp2matrixL1(1+(48*(i-1)):(48*i),5:7);
        obsRespErr = Exp2matrixL1(1+(48*(i-1)):(48*i),12);
    else
        newTestProbe = Exp2matrixL2(1+(48*(i-1)):(48*i),5:7);
        obsRespErr = Exp2matrixL2(1+(48*(i-1)):(48*i),12);
    end
    optTwo(i,:) = run(gs,problem);
end

% Median and mean parameters for full reference directions
med_optTwo_params(1,1) = median(optTwo(:,1));
med_optTwo_params(1,2) = median(optTwo(:,2));
mean_optTwo_params(1,1) = mean(optTwo(:,1));
mean_optTwo_params(1,2) = mean(optTwo(:,2));


refDirectionFull = 0;

for i=17:32
    if layout == 1
        newTestProbe = Exp2matrixL1(1+(48*(i-1)):(48*i),5:7);
        obsRespErr = Exp2matrixL1(1+(48*(i-1)):(48*i),12);
    else
        newTestProbe = Exp2matrixL2(1+(48*(i-1)):(48*i),5:7);
        obsRespErr = Exp2matrixL2(1+(48*(i-1)):(48*i),12);
    end
    optOne(i-16,:) = run(gs,problem);
end

% Median and mean parameters for 0 and 180 reference direction
med_optOne_params(1,1) = median(optOne(:,1));
med_optOne_params(1,2) = median(optOne(:,2));
mean_optOne_params(1,1) = mean(optOne(:,1));
mean_optOne_params(1,2) = mean(optOne(:,2));

for i=1:10000
    indx = randi([1 16]);
    bs_optTwo(i,:) = optTwo(indx,:);
    bs_optOne(i,:) = optOne(indx,:);
end

% Median and mean parameters for full reference directions (bootstrap)
bsmed_optTwo_params{p}(1,1) = median(bs_optTwo(:,1));
bsmed_optTwo_params{p}(1,2) = median(bs_optTwo(:,2));
bsmean_optTwo_params{p}(1,1) = mean(bs_optTwo(:,1));
bsmean_optTwo_params{p}(1,2) = mean(bs_optTwo(:,2));

% Median and mean parameters for 0 and 180 reference direction (bootstrap)
bsmed_optOne_params{p}(1,1) = median(bs_optOne(:,1));
bsmed_optOne_params{p}(1,2) = median(bs_optOne(:,2));
bsmean_optOne_params{p}(1,1) = mean(bs_optOne(:,1));
bsmean_optOne_params{p}(1,2) = mean(bs_optOne(:,2));

end
